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Abstract 

Using the Monte Carlo simulation method for bosonic reaction-diffusion systems introduced re- 
cently [S.-C. Park, Phys. Rev. E 72, 036111 (2005)], one dimensional bosonic models are studied 
and compared to the corresponding Langevin equations derived from the coherent state path in- 
tegral formalism. For the single species annihilation model, the exact asymptotic form of the 
correlation functions is conjectured and the full equivalence of the (discrete variable) master equa- 
tion and the (continuous variable) Langevin equation is confirmed numerically. We also investigate 
the cyclically coupled model of bosons which is related to the pair contact process with diffusion 
(PCPD). From the path integral formalism, Langevin equations which are expected to describe 
the critical behavior of the PCPD are derived and compared to the Monte Carlo simulations of the 
discrete model. 
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I. INTRODUCTION 



The reaction-diffusion (RD) systems have played a paradigmatic role in studying certain 
physical, chemical, and biological systems [lj. In the study of the RD systems on a lattice via 
Monte Carlo (MC) simulations, particles are usually assigned hard core exclusion property. 
On the other hand, the renormalization group (RG) calculations which have been successfully 
applied to several RD systems are often performed with boson systems 

sag 

. Hence, the 

comparison of the numerical studies with the RG calculations can sometimes become a 
nontrivial issue. 

There are two ways to bridge this gap between numerical and analytical studies. One 
is to make a path integral formula for hard core particles which is suitable for the RG 
calculations. This path has indeed been sought and some formalisms are suggested j^O,!?!. 
The other is to find a numerical method that would simulate boson systems. In this context, 
numerical integration studies of equivalent Langevin equations to the boson systems have 
been performed, toojo, . However, it is not always possible to find an equivalent 

Langevin equation [l2| and hence the applicability of this approach is somewhat restricted. 
Therefore, another numerical method is called for. 

Recently, a genera! algorithm to simulate the bosonic RD systems was proposed Q. 
Section |H] is devoted to a heuristic explanation of this algorithm to simulate general bosonic 
RD systems. In Sec. IIII1 the numerical method is applied to two bosonic RD systems. First, 
the single species annihilation model is studied with the emphasis on the pair correlation 
functions. We conjecture an exact asymptotic behavior of these quantities. We then present 
the numerical comparison of the discrete model to Langevin equation of continuous variables. 
Then, the cyclically coupled model of bosons is introduced and Langevin equations for this 
model with/without bias are derived from the well-trodden path integral formalism and 
compared to MC simulations. Section HVl summarizes the work. 



II. ALGORITHM 
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that is suitable for MC simulations 



This section explains the method proposed in Ref. 
of bosonic RD systems. After describing how single species boson systems can be simulated, 
a brief remark regarding the generalization to multiple species will be followed. 
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The reaction dynamics of diffusing bosons is represented as 

nA^(n + m)A, (f) 

where n > 0, m > —n (m / 0), and X nm is the transition rate. Each particle diffuses with 
rate D on a. d dimensional hypercubic lattice. The periodic boundary conditions are always 
assumed, but other boundary conditions do not limit the validity of the algorithm below. 
Configurations are specified by the occupation number p x (> 0) at each lattice point x. A 
configuration is denoted as {p} which means {p x |x £ L d }, where L d stands for the set of 
lattice points and the cardinality of L d is L d . 

The master equation which describes stochastic processes modeled by Eq. ([TJ takes the 
form 

BP / \ 

— = D ((Px + l)E x ,y ~ Px) P 

(x,y) ^ (2) 

+ ^ ^ nm {/(Px ~ m ' n )^x,m - /(Px, n) \ P, 

n,m x 

where P = P({p},t) is the probability with which the configuration of the system is {p} 
at time t, (x, y) means the nearest neighbor pair (x, y e L d ), /(p x , n) = (p x !)/(Px — n)\ is 
the number of ordered n-tuples at site x of the configuration {p}, and P Xjy and Cx, m are 
operators affecting P({p},t) such that 

P x , y P = P({--- ,p x + l,p y -l, ••■};£), 
C^P = P({- • • ,p x -m, ••■};*). 

The master equation implies that the average number of transition events for the config- 
uration {p} during infinitesimal time interval dt is 



E(dt,{ P }) = dt ( 2dP5n,l + Yl X ^ m J /(Px, 

x,n \ m / 

= dt^2 I 2d D5 nA + ^n!A nm J g-(p x , 

x.n \ m / 



(4) 

n), 



where g(p x ,n) = /(p x , n)/n! = is the number of (nonordered) n-tuples at site x. The 
first line of Eq. (J3J follows the usual convention in the field theoretical study of boson 
systems and the second line is introduced to save memories in actual simulations. For a 
later purpose, we introduce a model dependent function fo(p x ,n) = e n #(p x , n), where e n 
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takes 1 (0) if D5 Ut i + J2 m ^n,m is nonzero (zero). The meaning of e n is straightforward; we 
have only to consider the dynamics with nonzero transition rate. 

The algorithm starts by selecting one of n-tuples at any site, randomly. The simplest 
way to implement the selection is as follows: First a site x is picked up with probability 
JV X /M, where iV x = ^2 n h(p x ,n) is the number of accessible states (NAS) at site x and 
M = ^2 X N X is the total number of accessible states (TNAS). Then, n is chosen with 
probability h(p x , n)/N x . In this procedure, the array of the number of particles at all sites, 
say p[ ] (p[x] = p x ), is necessary. However, it is not efficient as there are too many floating 
number calculations. For a faster performance we introduce two more arrays, say list[ 
and act[ ][ ]. The array list[ ] refers to the location of any n-tuple. Each element of list[ 
takes the form (x, £), where x is a site index and £ lies between 1 and the NAS at site x. 
From £ and the array p[ ], which n-tuple is referred to by the array list[ ] is determined. If 
£ < h(p x , 0), then n = is implied. Else if I < h(p x , 0) + h(p x , 1), n — 1 is meant. Else if 
£ < h(p X: 0) + h(p x , 1) + h(p x , 2), £ indicates one of pairs at site x, and so on. In case the 
TNAS in the system is M, the size of list [ ] is M and all elements of list[ ] should satisfy that 
list [p] 7^ list [5] if p 7^ q (1 < p, q < M). Hence, the random selection of an integer between 
1 and M is equivalent to choosing one of all n-tuples with an equal probability. The array 
act is the inverse of the list. In other words, list [s] = (x, £) corresponds to act[x][£] = s. It is 
clear that these two selecting mechanisms are equivalent in the statistical sense. 

After choosing x and n, the reaction nA — * (n + m)A occurs with probability n\\ nm At for 
all m, where At is a configuration independent time difference. Provided n — 1 is selected, 
a particle at x hops to one of the nearest neighbors with probability DAt. To make the 
transition probability meaningful, At is made to satisfy 



for all n. After this update, time increases by At/M. On average, this algorithm generates 
E(At, {p}) transition events during At. 

For systems with k species, all we have to do is to modify the NAS at site x in such a 
way that 





k 





i=l n 



ni,...,n k 
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where the first (second) terms are from the dynamics in which n particles of same species (jij 
particles of each jth species) are involved. For instance, for the pair annihilation of different 
species, so-called A + B — > reaction, the second term of Eq. (JBJ) becomes P a,xPb,x- Except 
this modification, all other steps are the same as in the single species case. 

Equipped with the numerical methods, Sec. 11111 studies some bosonic RD systems which 
show scaling behavior. 



III. APPLICATIONS 



A. single species annihilation model 

The first example is the one dimensional single species annihilation model which corre- 
sponds to X nm = unless n = 2 and m = —2. For convenience, we set D = | and X2-2 — X. 
The decaying behavior of the particle density was studied in Ref . . This section studies 
the correlation function M(r; t) which is defined as 

1 L 

lim - J2(p x (t)p x+r (t)) if r ^ 0, 



M(r;t) 



L^oo L 

" L 



(7) 



lim ^(A^XA^)-!)) ifr = 0, 

x=l 

where (. . .) means the average over all independent realizations. Using the boson operators 
in Ref. 0,0,0], M(r;t) can be rewritten as ^ J2 x ( a x a x+r) ■ 

The correlation functions for the annihilation model of hard core particles with annihi- 
lation probability p were studied in Ref. jig . The asymptotic behavior of the correlation 
function is conjectured as [3] 

^^(iS^H^)' (8) 

with c = 3.4 ± 0.2. Note that M r (t) is not to be confused with M(r;t); M r {t) and M(r;t) 
are defined in the hard core and boson models, respectively. 

In fact, the exact value of c can be deduced from the differential equation 

^ = -2p Ml (t), (9) 

which relates the time derivative of the density p(t) to the correlation function with r = 1. 
Since p(t) ~ for any finite p in the asymptotic regime, it is easy to deduce that 



c = 7r if the asymptotic behavior of the correlation function takes the form of Eq. (JHJ. This 
value is compatible with the numerical estimation in Ref. Q|. 

By the same token, we can conjecture how the correlation function M(r; t) behaves asymp- 
totically from the equation 

dp(t) - -2AM(0;f). (10) 



dt 

If M(r; t) takes the similar form to Eq. (JHJ and since p(t) decays as 1/ y/Airt in the asymptotic 
regime for any nonzero value of A [l3|], one can deduce 

1 



M(r-t) ~ M^t) = 



7T 



r + 



(11) 



(4?rt) 3 / 2 V ' A, 

for all r > 0. As far as we are aware of, the correlation functions of the boson annihilation 
model have not been studied before. If D ^ |, the correlation function can be found by 
changing t ^ 2Dt and Ah>A/ (2-D). Since the boson model with infinite A is equivalent to 
the hard core particle model with p = 1 which is exactly soluble, Eq. (fTTj) becomes exact in 
this limit; see Eq. (jHJ. 

In the following, we will check the validity of Eq. (fTTj) for finite A and nonzero r via 
MC simulations. Initially, particles are distributed according to the uncorrelated Poisson 
distribution with average density po — 1. During simulations, we measured M(r;t) for 
r = 2°, 2 2 , 2 4 , and 2 6 up to t = 10 5 . The system size is 2 18 and around 2.5 x 10 5 independent 
samples are collected for both cases of A = 1 and ~. Figure [T] shows that M(r;t) takes the 
conjectured asymptotic form (jTTJ. 
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FIG. 1: Semi-log plots of M(r; t)/-Mas(r; t) as a function of t for (a) A = 1 and (b) A 
curves converge to 1 as t goes to infinity. 



i. All 



6 



The MC simulations of bosonic RD systems can confirm the equivalence between the 
(discrete) microscopic models and (continuous) Langevin equations, if exists. From the 
coherent state path integral representation of the bosonic systems |2[, Langevin equation 
can be derived in case each reaction involves at most two particles. Since the reaction 
of boson annihilation model requires two particles, one can write down Langevin equation 
which reads (ltd interpretation is employed) 

da x = dt(DV 2 x a x - 2Xa 2 x ) + iV2Xa x dW x , (12) 

where a x is a complex stochastic random variable whose average is the mean number of 
particles at site x, V 2 X is the lattice Laplacian defined as V 2 x f(x) = f(x+l)+f(x — 1) — 2f(x), 
i is the imaginary number, and W x is a Wiener process with (dW x dW x /) = dt5 x ^ x i. Initially, 
a x takes the value of p which is the initial density of the uncorrelated Poisson distribution 
used in the MC simulation. 

This equation is integrated using Euler scheme with At = 2.5 x 10~ 5 and the system size 
of 2 15 . In Fig. |21 numerical integration results for A = | are shown with comparison to MC 
simulations. Within statistical error, these two approaches yield the same results. Since the 
deviation from the mean field solution is evident, Langevin equation in the observation time 
properly appreciates the effect of noise. Hence, we believe that Fig. |21 shows the equivalence 
of two approaches for the annihilation model. Needless to say, the numerical integration of 
Langevin equation is a much harder job than the Monte Carlo simulation. 




1 2 3 4 5 

t 

FIG. 2: Plots of p(t) obtained from MC simulations (lines) and numerical integrations of Langevin 
equation (symbols) starting from the initial density pq. The broken line without symbols is the 
mean field solution of Eq. (|12|) . 
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B. cyclically coupled model 



The absorbing phase transition has been extensively studied as a prototype of the nonequi- 
librium critical phenomena Q|. The RG based on the boson systems has been applied 
successfully especially to the directed percolation (DP) universality class. Recently, the 
particle number probability distribution for boson models belonging to the DP class was 
studied numerically and the RG prediction was confirmed again |17| . 

On the other hand, the pair contact process with diffusion (PCPD) defies any numerical 
and analytical conclusions to date [3]. Although the driven PCPD (DPCPD) studied in 
Ref. Q| seems to conclude that the PCPD forms a different universality class from the DP, 
recent extensive numerical study 2(| revives the scenario that the PCPD will eventually be 
found to belong to the DP class with a huge corrections to scaling. Still, the universality 
classification for the one dimensional PCPD is yet to be settled unambiguously. 

To make matters worse, the recent RG study shows that the field theory starting from the 
single species master equation is not viable |2l|, which was also anticipated independently in 
Ref. 19] . As both works conclude, the field theory should account for the multispecies nature 
of the PCPD properly. Following this instruction, multi component Langevin equations with 



real random variables are introduced and studied in Ref. [22[ to find a viable field theory 
for the PCPD. We will take a slightly different path and ask whether we can find a viable 
field theory for the PCPD, in this section. 

Since the PCPD involves two independent "excitations" such as particles and pairs, it is 
natural to generalize to a two species model which captures the main physics of the PCPD. 
This type of two species model with hard core particles was introduced and studied in Ref. 

. This section introduces a bosonic variant and studies it using both MC simulations 
and Langevin equations. 

The model which will be called the cyclically coupled (CC) model is defined as follows: 
There are two species, say A and B. Each species diffuses with rate Da and Db, respectively. 
Each B particle is annihilated (B — ■> 0) with rate 5, branches another B particle (B — > 2B) 
with rate a, and mutates into two A particles (B — > 2A) with rate fi. Every pair of B 
particles at the same site can be coagulated (2B — > B) with rate A. Every pair of A 
particles produces a B particle and is removed (2A — > B) with rate r. The A (B) particles 
have a connection, if not a exact mapping, to the isolated particles (pairs) in the PCPD. 



S 
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FIG. 3: Semi-log plot of R(t) vs t for the CC with the relative bias near criticality. At criticality 
(Pc = 0.3751), clear logarithmic behavior is observed as in Ref. Inset: A plot of R(t)/\n(t) vs 
t in the semi-log scales. 

Since the PCPD as well as the CC suffers from the strong corrections to scaling, it is 
nontrivial to show directly by MC simulations that the CC and the PCPD should belong 
to the same universality class. Fortunately, we have an alternative to check the equivalence 
of the CC and the PCPD in the sense of the universality. If the relative bias between 
two species in the CC in one dimension triggers the mean field scaling with logarithmic 
corrections as happens in the DPCPD 19], it is reasonable to conclude that the CC and the 
PCPD share the critical behavior. 

The transition events of the CC with a relative bias in one dimension are almost same 
as those of the CC above except that A particles hop only to the right with rate 1. For a 
numerical study, we set D B — 0.1, p, — 0.2, r = 0.5, 5 = 2A = 0.6 x p, and a = 0.6 x (1 — p) 
with a tuning parameter p. Since only A particles diffuse in a biased manner, the relative 
bias between different species can not be gauged away by the Galilean transformation. 
Figure El shows that R(t) (= A(t)/ B(t)) which is a ratio of two densities at time t behaves 
logarithmically at criticality. Combining with the observation that A(t) ~ t~ 05 at criticality 
with possible logarithmic corrections (not shown), the CC with the bias shows the same 
critical behavior as the DPCPD, which confirms the equivalence of the CC to the PCPD in 
the sense of the universality. Accordingly, Langevin equations which are equivalent to the 
CC are supposed to describe the critical behavior of the PCPD. 

Following standard path integral formalism [2], one can derive the action of the CC, which 
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2 4 6 8 10 6 8 10 6 8 10 

t t t 

FIG. 4: (a) The densities of each species for both the biased and unbiased CC as a function of t 
from the MC simulations (lines) and numerical integration (symbols) of Langevin equations (jl4[>. 
For comparison, mean field solutions are also shown, (b) Close-up of the interval 6 < t < 10 for 
A(t) in (a), (c) Close-up of the same interval as in (b), but the plots are for B(t). 

reads 

£ = d x [d t a x - D A V 2 x a x + vd\\a x - 2{ib x + 2ra x ] 

+ b x [d t b x - D B W 2 x b x - rb x + Xbl - ra x ) (13) 

- hl(2ab x - 2\b 2 x ) - l -al{2^b x - 2ra 2 x ), 

where the average of the field a x {b x ) corresponds to the density of species A (B) at site x 
and r = a — fi — 5. Along the parallel direction denoted as ||, A particles hop to the right 
(left) with rate Da + v/2 (Da — v/2). Since the number of barred fields does not exceed 
two in each term, one can write down the equivalent Langevin equations to the action (|13|). 
which read 

da x = dt(D A V 2 x a x - vd\\a x + 2[ib x - 2ra 2 x ) + - 2ra 2 x dW x , (14a) 

db x = dt(D B V 2 x b x + rb x - \b 2 x + ra 2 x ) + y/2cxb x - 2\b 2 x dV x , (14b) 

where W x and V x are independent Wiener processes. 

In Fig. El we compare the MC simulations of the CC with the numerical integrations of 
Langevin equations (|14j) at p = 0.29. Initially, a x and b x are set to 1. The system size for 
the numerical integration is 2 15 and around 50 samples are independently generated with 
At = 2.5 x 10~ 5 . Up to t = 10, the difference between the unbiased and biased cases is 
minute, but, within statistical errors, the behavior of two cases can be discerned from each 
other. In other words, we showed that Eqs. (fbH) are equivalent to the CC with/without 
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bias. Although we compared two approaches just for one set of parameter values, the full 
equivalence for all parameter values is still expected. 

In summary, we showed that the CC and the PCPD share the critical behavior. Then, 
we found Langevin equations which are equivalent to the CC. From these two observation, 
we can say that Langevin equations (|14|) with complex random variables a and b show the 
same critical behavior as the PCPD. 

Although we found the representative Langevin equations for the PCPD, it is not obvious 
whether these equations with naive continuum limit can serve as a properly coarse-grained 
field theory for the PCPD. Besides, we are not sure whether Eqs. (I14j) contain all relevant 
(or sometimes dangerously irrelevant) terms. For example, the reaction A + B — > which is 
absent in our model can be generated by a chain of reactions. It is of no difficulty to write 
down Langevin equations with the pair annihilation of different species. However, what will 
happen if we include the reaction 3A — > which prohibits writing down Langevin equations 
like Eqs. (I14J) ? If this reaction is also important in whichever sense (relevant or dangerously 
irrelevant), terms with only a and a in the action take exactly the same form as those in 
Ref. (2]] . Hence, it seems that the difficulty found in Ref. 21 1 still remains even in the 
multi component Langevin equations studied here. We only hope that this study can be a 
starting point of the field theoretical understanding of the PCPD in the future. 



IV. SUMMARY 

To summarize, using the algorithm proposed in Ref. and generalized one to the 
multispecies models, the single species annihilation and the cyclically coupled models are 
studied. 

For the single species annihilation model, the exact asymptotic form of the correlation 
functions is conjectured and numerically confirmed. In addition, the equivalence of Langevin 
equation derived from the coherent state path integral formalism to the discrete boson model 
is affirmed. From the cyclically coupled model of bosons, we derive Langevin equations for 
both biased and unbiased cases. By simulating discrete models and integrating the Langevin 
equations numerically, these continuum equations are indirectly shown to describe the critical 
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behavior of the PCPD and the DPCPD. 
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